Discrete element modeling of particles sphericity effect on sand direct shear performance

Particle surface morphology is an important factor influencing sand structure and mechanical properties. In this study, the effect of sand particle sphericity on sand direct shear performance is investigated by using the discrete element method (DEM). Two ways are adapted to simulate different approaching methods from round particles to irregular sand. The macroresponse shows that irregular sand has a higher shear strength at lower normal stress than round particles. The shape of the particle has less influence on shear strength at higher normal stress. The irregular shape of sand leads to an increase in the shear band proportion. However, the shear band proportion is not related to the sphericity. Under all conditions, particles within the shear band have a larger average rotation angle than those outside the shear band. When the particle shape approaches round (regardless of the round particle proportion and particle shape), the average rotation angle of particles within and without shear bands increase, while the coordinate number and contact anisotropy decrease.


Direct shear test modeling and parameters calibration.
To benchmark the results from the corresponding numerical simulation, direct shear tests are carried out in parallel to calibrate the PFC input parameters. Clean and uniform sand with a specific gravity of 2.65 was chosen as the basic material in this research. The particle size distribution of the soil for the experimental and numerical tests is listed in Fig. 4. As the shear box width should be more than 10 times the maximum particle size and the shear box height should be more than 6 times the maximum particle size 29 , the shear box's dimensions for both the experimental direct shear test and numerical shear test was 40 mm * 40 mm. Sand was poured into the shear box with a dry density of 1.59 g/cm 3 . Loading speeding was set as 1% strain per minute. Direct shear tests were performed under normal stresses of  www.nature.com/scientificreports/ 50 kPa, 100 kPa, and 200 kPa. All tests were conducted under the guidance of the ASTM D3080-04 and Chinese standard GB/T 50123-2013. For the numerical direct shear test, the shear box is modeled by the enclosed wall.
The top wall is treated as a servo wall with a vertical velocity until the normal stress reaches the target value of 50 kPa, 100 kPa, or 200 kPa. When shear velocity is conducted at different velocities in the PFC, shear result differences are indistinct when the shear velocity is less than 0.1 m/s. This velocity is much different from the experimental test, as particle contacts may dissipate kinetic energy due to the damping acting 8 . Thus, the lower part of the shear box is made up of three walls, which are given a constant speed of 0.005 m/s, while the upper part of the shear box remains stationary. The stress on the vertical walls is recorded as shear stress, which is similar to the experimental test. It should be noted during the numerical simulation of the shear test that a horizontal plate is added to avoid particle leakage. The schematic diagram of conducted test is illustrated in Fig. 5. The contact model between the cohesionless sand is applied by the rolling resistance linear model 30 . The rolling resistance model is a linear-based model that adds a rolling resistance mechanism. The model incorporates a torque acting on particle contact to counteract rolling motion and to dissipate energy during the relative rotation 31,32 . The rolling resistance model can offer a simple method to simulate particle shape-like behavior 33 .    www.nature.com/scientificreports/ Following the proposed method by Lu et al. 8 , Young's modulus and friction coefficient are roughly calibrated. Then, a numerical direct shear test is conducted to match the experimental data. Figure 6 presents the shear stress-strain curve of the experimental test and numerical simulation for clean sand samples. The input parameters of the shear test are listed in Table 1.

Results and discussion
Shear strength of different particle assembles. The shear strength of sandy soils is mainly influenced of friction and interlocking structures between particle surfaces. Even under the same test conditions, the shear stress-displacement curves also present different trends. Figure 7 presents the shear stress-displacement curve for all cases under 50 kPa, 100 kPa and 200 kPa normal stress. At a lower normal stress (50 kPa), round particles show the lowest shear stress when compared with mixed soil or with a typical shape clump. There is a clear increasing trend of shear stress for the mixed soil. However, this increasing trend is unrelated to the round particle proportion. The medium (50% round disc 50% irregular clump), subirregular (25% round disc 75%    www.nature.com/scientificreports/ irregular clump), and irregular (100% irregular clump) soil present similar shear stress curves, while subround (75% round disc 25% irregular clump) shows smaller stress. For the single assigned type of clump, the shear stress increases as the S value decreases 34 . Even though particles have similar grain sizes and gradations, irregular particle shapes lead to a higher shear strength. At higher normal stresses (100 kPa and 200 kPa), the shear strength of mixed soil is still much higher than that of round particles. Even a small proportion (75% round disc 25% irregular clump) of irregular soil can significantly increase the soil shear strength. However, for the single assigned type of clump, particles with a larger S value exhibit little improvement in shear strength at higher normal stress. This indicates that when the particle approaches a round shape, its shape has less influence on shear strength at higher normal stress.

Shear band behavior.
When shear advances, the movement of coarse-grained particles results in particle relocation. The stress redistribution is adjusted by changing the shear displacement. The macromechanical shear performance of sand is the result of the microstructure's evolution during shear 35 . Figure 8 presents the www.nature.com/scientificreports/ irregular clump particle displacement vector at 50 kPa with shear displacements of 1 mm, 5 mm, and 10 mm. Different colors represent different displacement values, while the vector indicates the particle movement path. Particles in the lower shear box show larger horizontal displacement because they move with the lower shear box. Particles in the upper shear box exhibit two directions of movement: particles are forced to move downward at the top, while particles near the shear plane move horizontally. Particles near the shear plane show different displacement behaviors, but all present a bell-shaped curve. Soil particles with large displacement near the shear plane can be named shear bands or shear zones 36 . The localization of shear deformations and shear band formation contribute to the destabilization of soil. However, the traditional direct shear test mainly focuses on the macroscopic mechanical response. It is almost impossible to directly observe the behavior of particles in the shear zone. Thus, we quantitatively calculated the horizontal displacement of all particles in the shear box. Figure 9 presents the distribution of the soil particles' horizontal displacement for all cases under 50 kPa normal stress with shear displacements of 1 mm, 5 mm and 10 mm. The horizontal axis shows the particle's horizontal displacement, while the vertical axis shows the particle's vertical position in the shear box. The closer the particles are to the shear plane, the larger the horizontal displacement. Particles in the lower shear box mobilize with the shear box simultaneously; thus, the displacement of most particles in the lower shear box is similar to their horizontal movement.
To visualize shear band shape, horizontal displacement (L i ) of the sand was recorded. If Li/Ls ≥ 0.1 (Ls is shear displacement at the moment), this soil particle can be regarded as being within the shear band. Details are listed as follows.
when soil particles locat in the upper shear box  www.nature.com/scientificreports/ Then, the position information of the shear band particle was used to reconstruct the shape of the shear band in the scatter plot. Figure 10 presents the reconstruction images of the shear band distribution for irregular sand.
The shear band is usually simplified to have a rectangular shape. Formation of a shear band is universal. However, the shear band in the shear box has a fusiform shape with a large area in the middle. Zhou et al. 37 linked shear band width with rotation. Samples expanded with the development and reached stability. Vangla and Latha 38 estimated the shear band thickness from the initial and final profiles of colored sand columns. Although the shear band can be visualized, a quantitative definition of the shear band range is still lacking. To quantitatively analyze the shear band performance, the shear band proportion defined as the ratio of shear band particles to the whole particles in the shear box. The calculated shear band proportion is listed in Fig. 11. The left figures present the shear band proportion of all cases from 5 to 10 mm shear displacement under 50 kPa, 100 kPa and 200 kPa. The dashed line and corresponding percentage show the shear band proportion, while different colored solid lines show each particle's assembled shear displacement. The right figures list details of the shear band proportion that have shear displacements of 5 mm and 10 mm. The shear band proportion is not a constant value for each specimen; it changes according to the shear displacement. The irregular shape of sand leads to an increase in the shear band proportion. When shear advances, the rotation and movement of particles gradually become stable. It can be likely concluded that the shear band area increases with increasing shear displacement at each normal stress. As the normal stress increases, more particles are compacted into a dense state, which hinders particle movement during shear. Thus, the shear band area also increases. The irregular particles have a larger shear band area than the round particles. However, the shear band proportion is not related to the sphericity. There is no clear change trend of the shear band proportion from round to irregular particles. It should be noted that a mixture containing irregular and round sand increases the shear proportion significantly, especially at higher normal stresses.
Particle rotation. Particle resistance to rotation is known to develop in particle contact and contribute to enhancing the shear strength. As the shear band position and outline have been quantitatively obtained in the previous section, we classify particles into three groups: particles in the shear band, particles outside the shear band, and overall particles in the shear box. The corresponding average rotation angles are obtained and calculated. Figure 12 presents the average particle rotation angle at shear displacements of 5 mm and 10 mm under 50 kPa, 100 kPa and 200 kPa normal stresses. In all conditions, regardless of shear displacement and normal stress, particles within the shear band have a larger average rotation angle than those outside the shear band. When the particle shape begins to round (regardless of the round particle proportion and particle shape), the average rotation angle of particles with and without shear bands increases. This indicates that an irregular particle shape reduces particle rotation. In addition, compared with outside shear band particles, the average particle rotation angle within the shear band presents larger differences for each case, indicating that the particle shape has a larger influence on shear band particles. When shear advances from a 5 mm to 10 mm shear displacement, more particles are disturbed. Thus, a larger rotation angle is obtained. For round particles, the average contact number is small, as only one contact can be generated between two individual particles. The interparticle structure has difficulty resisting the moment caused by the tangential contact force, leading to a higher particle rotation. For irregular soil, the average contact number between particles increases, and the interparticle structure is more stable. Thus, the moment caused by the tangential contact force can be transmitted between the particles without causing excessive particle rotation. The excessive rotation of particles leads to a lower shear strength, which is confirmed in Fig. 7. The interlocking structure of irregular sand surfaces restrains and reduces particle rotation, resulting in a higher load resistance. When considering the effect of normal stress on particle rotation, it can be concluded that a higher normal stress increases the overall particle rotation angle. However, in the shear band, particles rotate less at a higher normal stress. When the normal stress increases, particle movement needs to overcome more external resistance, which leads to the reduction of particle rotation in the shear band. www.nature.com/scientificreports/ Particle contact behavior. In the PFC simulation, particles interact with others at contacts by generating internal forces. Individual particles are treated as rigid bodies in PFC. Thus, all deformations occur at the contacts. Though granular materials can be regarded as interacting particles, the particle contact force and orientation can determine the microstructure of the soil 39 . For the sand particle, forces are conveyed to one another by their contacts. Figure 13 presents the contact force chain distribution for round particles and irregular clumps before and after the test under 50 kPa normal stress. Force chain thicknesses are proportional to their magnitude. Before shearing, samples are subjected to compression force, as all samples are contained by the shear box. Because the shear box is subjected to vertical force, the compression force in the perpendicular direction 18  www.nature.com/scientificreports/ is higher than that in the horizontal direction. When the lower box starts to move, sand particles near the moving wall are forced to move. The contact force redistributes to be more diagonally oriented from the lower lefthand corner to the upper right-hand corner. In addition, the contact force increases with higher angularity 40   www.nature.com/scientificreports/ When particles approach irregularity, they are more discrete because more contacts can be formed between the irregular particles. Thus, the increase in contact force can be formed in irregular sand with a high interlocking structure, which is in line with Yang et al. 41 .
For the particles, the macrobehavior of shear strength mainly depends on the initial density and the normal stress, while the average contact number is another factor from a micro perspective 42 . The coordinate number is the number of particles in contact with the surrounding solids, and the average coordination number is the basic microscopic index of the granular material. The average coordination number (Z) is defined as where N c is particle contact number, N p is particles number. Figure 14 lists the particle coordinate numbers under 50 kPa and 200 kPa normal stresses. The coordinate number first increases and then quickly reaches a stable state. The coordinate number increases as the particle shape approaches an irregular shape. A good correlation can be found between the shear behavior and coordinate number. A higher coordinate number means that more contacts and interlocking structures are constructed, resulting in a steadier state. Thus, a larger shear strength can be obtained in more irregular sand 16,43 . Particles in the shear band have a larger coordinate number than the outside shear band particles. A higher normal stress also leads to a significant increase in the coordinate number. A higher coordination number leads to more sufficient contact between particles, thus leading to a particle skeleton system that can resist a higher external force.
The coordinate number shows the particle contact number performance and intergranular forces are also related to the contact orientation, which can provide information on the load transfer direction. It is known that the contact orientation of sand particles is closely related to its mechanical properties 44,45 . Under the shear test, the distribution of internal contact force depends on the contact orientation 46 . Thus, to obtain the sand particle orientation anisotropy performance during shear, the contact orientation distributions of the specimens before Figure 13. Contact force chain distribution before and after shear (A: round particle; B: irregular clump). www.nature.com/scientificreports/ and after shearing under 50 kPa normal stress are plotted in Fig. 15. At lower normal stress, round particles have an approximately isotropic contact direction before shear, but the contact direction becomes anisotropic when the particle shape becomes irregular. They present a windmill form with maximum contact numbers in the 0°, 90°, 180° and 270° directions. Because the specimen is confined by the rectangular shear box, only vertical stress is applied on the top plate of the shear box. The contact forces in the perpendicular and horizontal directions are dominant. In addition, the contact numbers of irregular sand particles are much higher than those of round sand. Because there is only one contact between two individual round particles, more contacts can be generated in the clumps due to their irregular shape. When the specimen is subjected to horizontal shear, the orientation of the contact force particles in the shear box begins to change. Because particles in the lower shear box are forced to move by the left plate, they are subjected to both vertical and horizontal forces simultaneously, leading to an increase in the contact number from 30° to 70° and 210° to 240°, respectively. This is consistent with the force chain distribution in Fig. 13, where the force chain presents a diagonal distribution after shear. Soil structure is the combination of soil composition, particle spatial arrangement and interparticle forces which results in the anisotropy 47 . The alignment direction of the particle leads to different contact forces and strengths. To quantitatively analyze the contact direction and the structural anisotropy of the corresponding material, a fabric tensor describing the particle arrangement and anisotropy is given as follows: where N is the contact number and n i and n j are the components of the unit branch vector in the i and j twodimensional directions, respectively.
In the two-dimensional plane, the fabric tensor can be expressed as: where cos θ (k) and sin θ (k) are k particle contact normal vectors in the i and j directions, respectively. F ij is a symmetric second-order tensor. According to the material plane stress state analysis method, the major fabric F 1 and minor fabric F 3 of F ij can be expressed as: Then, combining Eqs. (3) and (4), the following can be obtained: It is an anisotropy amplitude parameter that can be calculated to evaluate contact anisotropy. Figure 16 shows the contact anisotropy of each simulation at a 50 kPa normal stress. When particles are subjected to shear, the contact anisotropy value has a stable downward trend, highlighting the evolution of contact anisotropy. It can be seen from the left figure that round particles have the smallest contact anisotropy value, while irregular clumps have the maximum value. For the single shape particles, they all present anisotropic behavior. Contact anisotropy may not be directly linked to the S value as they (S value small, middle and large) all have similar values of contact anisotropy. On the other hand, the contact anisotropy value of the mixture containing round particles and irregular clumps in the right figure also shows a stable downward trend. However, with increasing the content of round particles, the contact anisotropy value decreases.

Conclusion
A discrete element simulation on direct shear performances of different particle shapes is presented. The sand particle shape is first captured by a microscope and then modeled by the particle flow code. Two methods are adapted to simulate different approaching methods for round to irregular sand. To analyze the micro behavior of sand under shear, we define the shear band proportion. This allows us to analyze the particle shear performance within or without a shear band. The shear stress-displacement relationship, particle displacement, shear band proportion, particle rotation, particle contact force chain, contact coordinate, and anisotropy are analyzed to evaluate the shear performance of sand with different shapes.
At lower normal stress, irregular sand shows higher shear strength than round particles. When the particle shape approaches a round figure, its shape has less influence on the shear strength at higher normal stress. The visualization of the shear band indicates that the shear band in the shear box has a fusiform shape with a thick shear band in the middle. The quantitative analysis of the shear band area proportion shows that the shear band proportion is not a constant value for each specimen; it changes according to the shear displacement. The irregular shape of sand leads to an increase in the shear band proportion. The mixture of irregular and round sand significantly increases the shear proportion, especially at higher normal stresses. Regardless of shear displacement and normal stress, particles within the shear band have a larger average rotation angle than those outside the shear band. Irregular particle shape reduces particle rotation. The coordinate number increases as  www.nature.com/scientificreports/ www.nature.com/scientificreports/ the particle shape forms into an irregular shape. A good correlation can be found between the shear behavior and coordinate number. A higher coordinate number means that more contacts and interlocking structures are constructed, resulting in a steadier state. Round particles have an approximately isotropic contact direction, but it changes to anisotropy when the particle shape gradually changes. www.nature.com/scientificreports/